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We perform fully relativistic calculations of binary neutron stars in corotating, circular orbit. While 
Newtonian gravity allows for a strict equilibrium, a relativistic binary system emits gravitational 
radiation, causing the system to lose energy and slowly spiral inwards. However, since inspiral 
occurs on a time scale much longer than the orbital period, we can treat the binary to be in 
quasiequilibrium. In this approximation, we integrate a subset of the Einstein equations coupled 
to the relativistic equation of hydrostatic equilibrium to solve the initial value problem for binaries 
of arbitrary separation. We adopt a polytropic equation of state to determine the structure and 
maximum mass of neutron stars in close binaries for polytropic indices n = 1, 1.5 and 2. We 
construct sequences of constant rest-mass and locate turning points along energy equilibrium curves 
to identify the onset of orbital instability. In particular, we locate the innermost stable circular 
orbit (ISCO) and its angular velocity. We construct the first contact binary systems in full general 
relativity. These arise whenever the equation of state is sufficiently soft (n ~ 1.5). A radial stability 
analysis reveals no tendency for neutron stars in close binaries to collapse to black holes prior to 
merger. 

PACS numbers: 04.20.Ex, 04.25.Dm, 04.30.Db, 04.40.Dg, 97.60.Jd 



I. INTRODUCTION 

Neutron star binaries are interesting for numerous rea- 
sons. Several neutron star binary systems are known to 
exist even within our own galaxy ||l| . For some of these 
systems (including PSR B1913-M6, B1534-hl2) general 
relativistic effects in the binary orbit have been measured 
to high precision |^,^ . Binary neutron stars are believed 
to be among the most promising sources of gravitational 
waves for detectors like LIGO, VIRGO and GEO. This 
circumstance has triggered multiple efforts to predict the 
gravitational waveform emitted during the inspiral and 
the final plunge of the two stars. More fundamentally, 
the two-body problem is one of the outstanding unsolved 
problems in classical general relativity. 

Considerable effort has gone into understanding bi- 
nary neutron stars. Most of this work has been per- 
formed within the framework of Newtonian hydrodynam- 
ics. Hachisu and Eriguchi Q constructed hydrostatic 
equilibrium of binaries in synchronized circular orbits. 
Rasio and Shapiro studied binary equilibrium con- 
figurations and their dynamical evolution, including the 
merger of the two stars. The coalescence of neutron star 
binaries has also been investigated by Shibata, Nakamura 



and Oohara Zhuge, Centrella and McMillan j?] and 
Ruffert, Janka and Schafer and other investigators. 

Many investigators have also studied the binary prob- 
lem within a post-Newtonian framework. As long as the 
stars are well separated they can be approximated by 
point sources. In this case hydrodynamical effects are 
neglected and the gravitational waveform can be calcu- 
lated to second post-Newtonian order (see ||] and ref- 
erences therein). Post- Newtonian calculations that do 
take into account hydrodynamical effects are also un- 
der way: Shibata and Taniguchi and Shibata |ll| ] 
have constructed equilibrium configurations and Oohara 
and Nakamura ||l2(| have studied binary coalescence. 
Lai Lai and Wiseman [Q and Lombardi, Rasio and 
Shapiro [Tsl] have constructed binary equilibrium config- 
urations in an ellipsoidal approximation. 

Fully general relativistic treatments of the problem 
are complicated by several factors, including the non- 
linearity of the partial differential equations and the re- 
quirement of very large computational resources to solve 
the coupled system. These simulations are currently 
only in their infanc y fl^ . Recently, Wilson, Math- 
ews and Marronetti [|I6[ (hereafter WMM) reported re- 
sults obtained with a relativistic numerical code. Their 
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code assumed several simplifying physical and mathe- 
matical approximations. Their results suggest that the 
central densities of the stars increase as the stars ap- 
proach each other and that massive neutron stars indi- 
vidually collapse to black holes prior to merger. WMM 
therefore find that in general relativity, the presence of 
a companion star and its tidal field tend to destabilize 
the stars in a binary system. This conclusion is op- 
posite to what is exp ected from Newtonian [jl^, post- 
Newtonian |13-l^,|l^, perturbative [|l9| and matched 
asymptotic expansion pO| , pT[ treatments of the problem. 
WMM also find that just prior to plunge and merger, 
their binary system has a total angular momentum too 
large to form a Kerr black hole (see the discussion in [p2t). 

In this paper we construct fully relativistic binary neu- 
tron stars in quasiequilibrium circular orbit ( "quasi" - 
equilibrium because these binaries are not strictly sta- 
tionary: because of the slow emission of gravitational 
radiation, general relativistic binaries cannot be in strict 
equilibrium). These models are interesting on their own 
right and provide initial data for future dynamical evolu- 
tion calculations. We study the structure of the neutron 
stars in these close binary systems and determine, for 
example, their maximum allowed equilibrium mass. In 
addition, we build quasiequilibrium binary sequences of 
constant rest-mass. These sequences approximate evolu- 
tionary trajectories of neutron star binaries undergoing 
slow inspiral via the generation of gravitational radiation. 
By locating the turning points in their total energy ver- 
sus separation curves, we can identify the onset of orbital 
instability at the innermost stable circular orbit (ISCO) 
and the orbital parameters at that critical radius. We 
have presented preliminary results in p^ , and analyzed 
the stability of these binaries in |24) . We do not find any 
evidence for a destabilization of neutron stars in close 
binaries. 

The purpose of this paper is to discuss details of our 
approximations, equations, and numerical method, and 
to present more complete results. The paper is organized 
as follows: in Section || we discuss all the underlying as- 
sumptions and approximations made in our calculations. 
In Section [I] we derive all the equations describing the 



quasiequilibrium of relativistic binary neutron stars. The 
numerical implementation of these equations is described 
in Section IV. We present results for several different 
polytropic equations of state in Section ^ and briefly 
summarize our findings in Section^. We also include an 
Appendix with tabulated data for some of our sequences. 



II. BASIC ASSUMPTIONS AND 
APPROXIMATIONS 

Throughout this paper we will assume that the two 
neutron stars have equal mass, are corotating in a circular 
orbit and that the matter obeys a polytropic equation of 
state. 



Choosing a polytropic equation of state permits a wide 
survey of models as a function of the stiffness of the 
equation of state and also simplifies the integration of 
the matter equation (|33|). However, polytropic equations 
could be easily replaced by more realistic cold equations 
of state. 

Restricting our analysis to stars with equal masses al- 
lows us to exploit spatial symmetry and solve the problem 
in just one octant in our Cartesian grid (see Section [II A 
below). However, generalizing our method to stars of 
unequal mass is straightforward. Nevertheless, it is in- 
teresting to note that all well determined masses of neu- 
tron stars in close binary systems have masses remarkably 
close to IAMq (see, for example, 0). Focussing on stars 
with equal mass may therefore be physically reasonable 
as well as numerically convenient. 

Demanding that the stars be corotating is a much 
less realistic assumption. Even if the stars in a binary 
started out corotating at a large separation, maintain- 
ing this corotation during inspiral would require a larger 
viscosity than is possible in neutron stars |^,^. In- 
stead, it is more likely that the circulation of the stars 
is conserved during inspiral. However, our assumption of 
corotation greatly simplifies the solution of the problem 
(see Section [II B) and it is appropriate to tackle this 
simpler case first. Even in Newtonian theory, the con- 
struction of nonsynchronous binaries is difhcult because 
of the unknown velocity field; only in ellipsoidal models 
can one build nonsynchronous as easily as synchronous 
binaries ||l5|] . Constructing more realistic sequences of 
constant circulation requires a dynamical treatment, as 
one marches inward from one radius to the next using 
the full coupled set of field and hydrodynamic evolution 
equations to guarantee conservation of circulation. 

In Newtonian gravity, a strict equilibrium solution for 
two such stars in a synchronized circular orbit always 
exists, except for very stiff equations of state (with n ~ 
1.5) near contact Since this solution is stationary, 
the hydrodynamical equations for the matter reduce to 
a single Bernoulli in tegral , which greatly simplifies the 
problem (see Section [II C ). 

Because of the emission of gravitational waves, a bi- 
nary in general relativity cannot be in strict equilib- 
rium. However, up to the innermost stable circular or- 
bit (ISCO), the timescale for orbital decay by radiation 
will be much longer than the orbital period, so that the 
binary can be considered to be in "quasiequilibrium" . 
This allows us to neglect both gravitational waves and 
wave-induced deviations from a circular orbit to good 
approximation. A similar approximation is often used 
in stellar evolution calculations: there the relevant evo- 
lution timescales are the nuclear or Kelvin-Helmholtz 
timescales, while the stars maintain (quasi) hydrostatic 
equilibrium on a dynamical timescale. 

We attempt to minimize the gravitational wave con- 
tent by choosing the spatial metric to be conformally 
flat, as in WMM (see also [p^). As wiU be shown in 
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Section III A , the field equations then reduce to a set of 



coupled, quasi-linear elliptic equations for the lapse, the 
shift and the conformal factor. If we neglect small devi- 
ations from circular orbit, the fluid flow is again station- 
ary, and the hydrodynamical equations ag ain re duce to a 
relativistic Bernoulli integral (see Section III B ). Solving 
these equations yields a valid solution to the initial value 
(constraint) equations, and an approximate solution to 
the full Einstein equations for an inspiraling binary at 
any given moment, prior to plunge. 

This conformal approximation has been carefully 
tested in Ref. |^8| for a single rotating star in station- 
ary equilibrium, which is the simplest numerical example 
in relativity for which the equilibrium solution deviates 
from conformal flatness. In Ref. ||2^ it was shown that 
by assuming conformal flatness, the resulting deviations 
from the exact solution were typically much smaller than 
1%, even for highly relativistic stars. 



III. BASIC EQUATIONS 

A. Field Equations 

To construct a numerical model of a binary system we 
employ the ADM decomposition of Einstein's equations 
of general relativity ||2^ . The derivation of our adopted 
equations closely follows the derivation in ||2^ for rotating 
stars. 

We write the metric in the general form 

ds^ = -a'^dt^ + J^j{dx' - Lu'dt){dx^ - ujHt). (1) 

Throughout the paper Latin indices will run from 1 to 3, 
whereas Greek indices will run from to 3. We also set 
G = c = 1 . By definition of the extrinsic curvature Kij , 
the three-metric 7^ satisfies the dynamical equation 



(2) 



where denotes the covariant derivative associated with 
7y . This equation can be decomposed into its trace 

dtln-f^^^ ^ ~aK - DiUj\ (3) 
where 7 = det "fij and K — K''^, and its trace- free part 
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-DiUjj ~ DjUJi + --jijDkOJk- (4) 

In the following we will choose maximal slicing so that 

K ^0. (5) 

We expect the gravitational wave content of the space- 
time to be small (see Section ||) , and we now want to use 
this expectation to simplify the problem. Unfortunately, 



the physical fields cannot be cleanly separated into freely 
specifiable dynamical degrees of freedom and dependent 
quantities, which are determined by the constraint equa- 
tions. However, such an identification is possible with the 
help of a conformal decomposition . We can therefore 
attempt to minimize the gravitational wave content of 
the (physical) spacetime by removing the dynamical (or 
"wave") degrees of freedom from the conformal fields. 
This can be achieved by choosing the three-metric 7^ to 
be conformally flat, so that "f~^^^jij = fij, where is 
the flat space metric. We will later use Cartesian coor- 
dinates, for which /y becomes the Kronecker delta 6ij. 
Note that this choice can always be made to find initial 
data on one time slice without any approximation. Our 
approximation lies in assuming that the metric will re- 
main conformally flat for all times during the inspiral. 
Eq. (0) then reduces to Q 



2aK,j = -DiLOj - DjLj, + -^ijDkUJ '. (6) 



We now write the metric as 



(7) 



where 5* is the conformal factor. The later is determined 
by the Hamiltonian constraint 



where the source term p is defined by 
p = n n'^Taf). 



(8) 



(9) 



Here is the normal vector to a i = const slice and Tap 
is the stress-energy tensor. For the metric (0), the Ricci 
scalar R in (0) reduces to 



(10) 



where is the fiat space Laplacian associated with . 
Inserting this into (0) we find 
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(11) 



Here we have transformed K'^^ according to 

^u^^io^'U^ (12) 

which, from Eq. (p|), now satisfies 



/ . . . . 2 . ■ i. 



(13) 



Inserting this expression into the momentum constraint 
DjK'^ = Snf (14) 

yields 
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Here the source term is given by 



(15) 



(16) 



This equation can be simplified by writing the shift vector 
as a sum of a vector and a gradient [B2| 
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(17) 



Eq. ( [Tsl ) can then be replaced by the two equations 

V^G" = 2Vj lTi{a^-^)K'^ - 16™^'''/ (18) 

and 

V^B = V^G'. (19) 

Imposing the full set of dynamical equations for the 
evolution of Kij would be inconsistent with Eq. (^ and 
our approximation that 7^ remains conformally flat at 
all times. However, in addition to Eq. (||) we can always 
require that the maximal slicing condition be preserved 
dtK = 0. Taking the trace of the time evolution equation 
for Kij together with Eq. (|l^) then yields an equation for 
the lapse 



Here the source term S is defined by 



(20) 



(21) 



Eqs. (0), (^, (|T|) and ^ together with the matter 
equations (see next Section) form a system of coupled, 
nonlinear elliptic equations, which has to be solved itera- 
tively. The boundary conditions follow from asymptotic 
flatness. Following Bowen [Q, the exterior solution to 
the field equations can be expanded in terms of multi- 
pole moments. We adopt as outer boundary conditions 
the fall-off behavior of the lowest order non-vanishing 
multipole moments. Because of the symmetries of the 
problem it is possible to solve it in only one octant of 
a Cartesian grid. The resulting boundary conditions on 
the coordinate planes together with the outer boundary 
conditions are summarized in Table 1. 



B. Matter Equations 

As we have discussed in Section we neglect wave- 
induced deviations from a strictly periodic, circular orbit, 
and also assume the stars to be corotating. In Cartesian 
coordinates we can choose the equatorial plane to be the 
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TABLE I. Boundary conditions for the outer boundaries 
[r 00) and on the coordinate planes in Cartesian coordi- 
nates. The equatorial plane is taken to be the y = plane 
and the stars are taken to be aligned with the z-axis. 



2/ = plane, so that the fluid four velocity then takes the 
form 



(22) 



where f2 is the constant angular velocity. We introduce 
a vector 

r = (0,z,0,-a:), (23) 

in terms of which the four velocity can also be written 

u" =w*(an" + f7r-w")- (24) 

Define v to be the relative velocity between the matter 
and a normal observer 
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au 



(l-f2)l/2 

Then, from = —1, we find 

For a perfect fluid the stress energy tensor is 

T-^f^ = (po + + P)u'^u'^ + Pg'-l', 



(25) 



(26) 



(27) 



where po is the rest-mass density, pi is the internal energy 
density and P is the pressure. The source term p in 
Eq. (^) can then be written 



Po + Pi + P „ 

P = — ] 2 ^' 

1 — w-^ 

the momentum source / in Eq. ( |l6| ) becomes 

{po + p^ + P) (nc-to^) 



J = 



(28) 



(29) 
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and S in Eq. (]21|) is given by 



(PO + Pz + P) 



1 



3P. 



(30) 



In order to describe the matter close to equilibrium we 
will use two of our basic assumptions. Neglecting devi- 
ations from a strictly periodic circular orbit and taking 
the two stars to be corotating is equivalent to assuming 
that the fluid four velocity is proportional to a Killing 
vector 



d 

1- — . 

dt d(b 



(31) 



In this approximation, the matter equations can be inte- 
grated analytically, which yields the relativistic Bernoulli 
integral (see, e.g., [§3) 



— = const. 
h 



Here h is the enthalpy 
h = exp 



dP 



For a polytropic equation of state 

D 1 + 1/n 

KPO ' , 



(32) 



(33) 



(34) 



where n is the polytropic constant and n the polytropic 
index, the enthalpy becomes 



PO +P;+P 
PO 



(35) 



It is very useful to introduce a dimensionless ratio 



P 

9 = — , 
Po 

in terms of which we can express 



(36) 



(l + (l + n)g) jn^'-cu') 
a 1 — v"^ 



p + 2S = K-^q- [ i±^i±^(l + 2v') + 5q 



(41) 
(42) 
(43) 



Note that physical dimensions enter our problem only 
through the polytropic constant k in the equation of 
state (p^). It is therefore useful to nondimensionalize all 
equations and eliminate k from the problem. This means 
that given the polytropic index n, we can solve the equa- 
tions once and use the results for arbitrary n. Since k"/^ 
has units of length we can introduce dimensionless coor- 
dinates t = K~"'^^t, X = and the same for y and 
z. The derivative operator scales as Vi — K"/^Vi, and 
the extrinsic curvature as K^^ = . The angular 
velocity Vl transforms according to fj = We also 
rescale B = kT'^/'^B and = k^"/^^' Putting terms 
together we find the Hamiltonian constraint 



2 " 



1 - 



(44) 



the lapse equation 
8 

+27r5*V (^(l + (« + l)g)^j-^^ + 5g) , (45) 
and the momentum constraint equations 



-16.vl/V^ + (^+"^'^(»e-c.') 
1 — 



(46) 



and 



P 



(37) 
(38) 
(39) 



Note that in the Newtonian limit we have g <C 1. In- 
serting the last three expressions together with Eqs. ( |25| ) 
and (pq) into Eq. (B2|) we find 



1 



9 = 



l + C 



l + n Va(l -i;2)i/2 



(40) 



where we have written the constant in Eq. ( p2|) as l + C . 
Also, we use q to rewrite the source terms (P8|-p0|) as 



Here we have used 



a — '^a. 



(47) 



(48) 



Eqs. (|44|-[17D together with ( }40| ) form a set of seven equa- 
tions for the seven unknowns a, G% B and q. More 
specifically, we have to find a solution to six coupled, 
quasi-linear elliptic equation for the gravitational fields, 
together with one algebraic equation for the matter. K^^ 
and Lo"^ in the above expressions can be expressed in terms 
of the unknowns with the help of Eqs. ( |l3|) and (|l7|). 
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C. The Newtonian Limit 



x = x/rB y = y/rB z = z/rs, (55) 



In this Section we will briefly show that in the New- 
tonian limit the above equations approach the expected 
form. In particular we expect 



(49) 



where $ is the Newtonian potential. Also, in the Newto- 
nian limit $, C, w ^ 1, so that (Eoh becomes 



71+ 1 V 2 



J_(c-<p + ln^{x^ + z^^ 

n+1 \ 2 ^ 



(50) 



Here we have used oj^ — (absence of frame dragging in 
the Newtonian limit). This limit, by Eq. (|l^), implies 
K'^ = 0. With g < 1, Eq. (||) now reduces to 



5 „n 



Identifying 



2 



(51) 



(52) 



yields, to leading order, the Poisson equation 

V^* = 47rpo- (53) 
Eq. ( ^ ) gives the same limit. 

IV. NUMERICAL METHOD 
A. Constructing Quasiequilibrium Models 

Corotating, equal mass binaries in circular orbits form 
a two-parameter family (just like single, uniformly ro- 
tating stars). A particular configuration is uniquely de- 
termined by two independent parameters. For compu- 
tational purposes it is particularly convenient to choose 
these parameters to be the maximum density gmax and 
the relative separation of the stars . 

As mentioned in Section [II, we choose the stars to 
orbit in the y = plane and to be aligned with the z-axis. 
In this case the surface of one star will intersect the z- 
axis at two different places. We will label the intersection 
closer to the origin of the coordinate system ta and the 
one further out fs- The ratio 



ZA = ta/tb 



(54) 



then parametrizes the relative separation of the stars. We 
can construct an algorithm for solving the gravitational 
and matter equations by modifying the algorithm used by 
several authors for single rotating stars |^,^,^. Mak- 
ing this algorithm stable requires rescaling the coordi- 
nates with respect to so that 



which means that the outer edge of the matter will always 
be at = 1. We also rescale 



(56) 



Eqs. (^) and (p^-]47|) are left unchanged, except that the 
matter source terms in (44-46|) have to be multiplied by 
and Vi has to be replaced by V j . This rescaling then 
allows us to determine rs as well as the the angular veloc- 
ity Cl and the matter constant C via an iteration process 
that uses Qmax and za as the two input parameters. 

The iteration scheme starts with an initial guess for 
the rest density distribution. We chose the density pro- 
file of an isolated, spherical star, i.e. we integrate the 
Tolman-Oppenheimer-Volkoff equations for the central 
density ^max and rescale the profile such that it fits be- 
tween rA — ZA and = 1. For this matter distribution 
we can then find a solution to the field equations (|4^-|4^) 
using a full approximation storage multigrid scheme (see, 
e.g., H). 

Once a solution to the field equations has converged to 
an adequate accuracy on the finest level of the grid hier- 
archy, we evaluate Eq. (|o|) at three different locations to 
find new values for the constants fi, C and fs as well as 
a new density distribution. To do so we first search for 
the maximum density along the z-axis [ [40[ and call this 
location fc- We can then evaluate Eq. pO| ) at the three 
points f^, fs and fc 



(1 + (n + l)q) ( - ^^{nz - w^f 



1/2 



(57) 



where we have used x =^ y = lu^ — u'' — on the z- 
axis. Note that at ta and the density vanishes q — 0. 
This set, on first sight, looks like three equations for the 
two unknows Q and C. However, changing the scaling 
parameter tb will also change the gravitational fields, so 
that a and ^> will implicitly depend on f^. We determine 
how a and 4* scale from the Newtonian limit. Rescaling 
the Poisson equation shows that the Newtonian potential 
$ scales with fg. Eqs. ( ^ ) and ( p2| ) therefore suggest 
that a and "if should be rescaled according to 



1- = (* 



(58) 



Inserting these scale relations into (^) then yields three 
equations for the three constants Cl, C and fs, which 
can be solved iteratively. Once the constants have been 
determined the new matter distribution can be calculated 
using Eq. (|^). 

The iteration can then be continued by finding the new 
fields for the new matter distribution. At each step we 
calculate the residuals of Eqs. (|44|-^7|) and integrate these 
over the numerical grid. We typically stop the iteration 
when the sum of these six integrated residuals is smaller 
than about 1 % of the estimated truncation error on the 
finest grid. 
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Once an iteration has been completed, we can calculate 
several physical quantities that characterize the configu- 
ration. The total rest-mass Mq tot is 



Im Jm 



-gcfx, (59) 



where the subscript M denotes integration over the sup- 
port of the matter and \/—g = In nondimensional 
form wc can therefore write 



Afo.tot = K-"/'Mn.tot = 
The total mass-energy (ADM mass) is 



o,tot = I a^^^u'q'^dH. (60) 

M 



— (f \/''i>(fS, = — - ( V2*d3^, (61) 
27r /oo 27r 



Using the Hamiltonian constraint ( |4^ ) this can be rewrit- 
ten 

1 



Mtot = / ■^-^k.jk'^d'x 

167r ' 



M 



(62) 



or, in nondimensional form, 
Mtot = K-"/'Mtot 

M 



167r 



1 - W2 



- q <rx. (63) 



Eq. (^ is the actual form we use to evaluate Mtot- The 
angular momentum is aligned with the y-axis and can be 
defined as 

j^^^ ^bill I x^k'^'d^S, = ^ [ x^Vik'^'d^x (64) 

(see, e.g., [^). This is the total angular momentum 
contained in the spacetime and includes both the orbital 
and spin angular momentum of the stars. Using Vj-ftT*''' = 
'^^'^DiK^^ as well as the momentum constraint (111), this 



can be rewritten 
Jtot = 

JM 



^^°{zf -xj'-)d^x^ [ •i/^^f.jCfd^x, 
Jm 



(65) 

where we have also used definition (p3|). Finally, we can 
substitute ( ^ ) for and write the angular momentum 
in the nondimensional form 



Ju 
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1 + (1 + n)g 



f.jCini^ -iu^)d^x, (66) 



where we have rescaled ^ according to ^ = S,/^b- 



In the following we will denote half the total rest-mass, 
mass and angular momentum by Mq = Mg, tot/2, M = 
Aftot/2 and J = Jtot/2. In the limit of large separation, 
Mq and M approach the corresponding values of isolated 
stars. 

Performing numerical simulations in three dimensions 
requires large computational resources. We have there- 
fore implemented our algorithm in a parallel environment 
using the DAGH infrastructure and run it both on 
the SP2 cluster at the Cornell Theory Center and the 
Origin2000 at the National Center for Supercomputing 
Applications at the University of Illinois. We typically 
use grids of (64)'^ or (128)'^ gridpoints, and run the code 
in parallel on 8 processors. DAGH has been developed as 
part of the Binary Black Hole Grand Challenge Project 
and is a package of routines and computational structures 
that allows for a convenient implementation of parallel 
applications on grid hierarchies. 



B. Constructing Quasiequilibrium Sequences 

In addition to constructing individual quasiequilibrium 
configurations, we can also build quasiequilibrium se- 
quences of constant rest- mass Mq. As we will discuss 
in Section these sequences provide approximate evo- 
lutionary tracks of inspiraling neutron star binaries. 

Our quasiequilibrium configurations are parametrized 
by their relative separation za and maximum density 
<Zmax- We therefore have to find a path through this two- 
dimensional parameter space along which Mq is constant. 
This can be achieved in several different ways. For ex- 
ample, for each separation za one could vary Omax until 
a configuration of mass Mq has been found Here 
we found it easier to start with a small (and hence only 
mildly relativistic) gmax for each za, and then increment 
9max in small steps keeping za constant. The results can 
be tabulated, and the procedure repeated for a different 
ZA- Once sufficient data have been collected one can then 
interpolate to a chosen rest-mass Mq. Note that for each 
za = const sequence we adjusted the outer boundary so 
that the number of gridzones covering the stars is the 
same for all separations. 

We have performed several tests to check our code. In 
two different regimes the results can be compared with 
known solutions: for small masses and weak fields we 
recover the Newtonian limit, and for large separations 
we approach the Oppenheimer-Volkoff spherical solution 
for each star and its near-by field. We have also checked 

"Til 



the fully relativistic identity |42 



dMtot = ndJt, 



(67) 



which holds along constant rest-mass sequences. To eval- 
uate Eq. (^7|), we have to take numerical differences be- 
tween integrals of very similar magnitude, so that their 
relative error was much larger than that of the individual 
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FIG. 1. Rest-density contours in the equatorial plane for 
a neutron star binary close to the ISCO. Each star has a 
rest-mass of Mo — 0.169, corresponding to a compaction in 
isolation of {M/R)oo = 0.175. The contours show isosurfaces 
of the rest-density in decreasing factors of 0.556. 



FIG. 2. Rest-Mass Mo versus maximum density pc for sep- 
arations ZA ~ 0.3 (bottom solid line), 0.2, 0.1 and 0.0 (top 
line). The dashed line is the Oppenheimer-Volkoff result for 
a n = 1.0 polytrope. The insert is a blow-up of the region 
around the maximum mass. 



integrals. Nevertheless, we found that this identity is sat- 
isfied typically to ^ 10% (except close to turning points, 
where the error due to the differentiation dominates). We 
expect that the numerical data presented in this paper 
are typically accurate to within a few percent, and are 
confident that our code correctly predicts qualitative fea- 
tures, like, for example, changes in the maximum allowed 
mass. 



V. RESULTS 

A. Sequences for n = 1.0 

In this Section we discuss configurations and sequences 
with a poly tropic index n = 1, representing a fairly stiff 
equation of state. This is a particularly interesting ex- 
ample, since realistic neutron stars are expected to be 
governed by equations of state of similar stiffness. Re- 
sults for n = 1.5 and n = 2 will be presented in Sec- 
tion VB. Numerical values in geometrized units can be 



obtained from our nondimensional "barred" quantities 
by multiplying with ap prop riate powers of k, according 
to Eqs. (|g), (H) and (||) (for example M = k"'/'^M, 
J = J, and po — n^^po). 

In Fig. 1 we show the density profile in the equatorial 
plane of a binary neutron star. Here za = 0.175, and 
the stars are close to the ISCO (see below). Each star 
has a rest-mass of Mq = 0.169, corresponding to a com- 
paction in isolation of {M/R)oo ~ 0.175. The contours 



show isosurfaces of the rest-density in decreasing factors 
of 0.556. The maximum compaction of a stable n = 1.0 
polytrope in isolation is {M/R)oo = 0.216, correspond- 
ing to a maximum rest-mass Mq = 0.180 and a maximum 
mass M = 0.164. 

In Fig. 2 we plot the rest-mass Mq versus the maximum 
density pc = p^'^^ + pf^^ for several different separations 
between za = 0.3 (roughly two stellar radii apart) and 
= (touching). As — > 1, we expect these curves 
to approach the spherical Oppenheimer-Volkoff (OV) re- 
sult, which we included as the dashed line in Fig. 2. Note, 
however, that the exact OV curve is computed from a 
one-dimensional ordinary differential equation with very 
high accuracy, while the binary configurations have been 
calculated on very coarse, three-dimensional numerical 
grids. From convergence tests we know that we system- 
atically underestimate masses, and accordingly, for large 
separations, we find masses slightly smaller than the cor- 
responding OV masses. All graphs lie within less than 2% 
of the OV curve, showing that the presence of a compan- 
ion star has only very little influence on the mass-density 
relationship. 

As we decrease the separation, the mass supported by a 
given central density pc increases slightly. In particular. 



the maximum rest-mass increases from MR 



= 0.179 



for ZA = 0.3 to M^'^'^ = 0.182 for stars in contact. This 
trend clearly suggests that the maximum allowed mass 
of neutron stars in close binaries is slightly larger than 
in isolation. This increase is caused in part by the ro- 
tation of the stars and in part by the tidal fields. More 
specifically, we find that the increase of the maximum al- 
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FIG. 3. Binding energy and angular momentum as a func- 
tion of the angular velocity for several different values of Mo • 
The curves are labeled by the compaction [M/ R)ao of the 
stars in isolation at infinity, starting with 0.05 and increasing 
in steps of 0.0025 up to 0.2. The maximum compaction of 
a stable, isolated, non-rotating n — 1.0 polytrope is 0.217. 
The upper label gives the orbital frequency for stars with a 
rest-mass of 1.5M0 



lowed mass is comparable to the corresponding increase 
of an isolated neutron star rotating with the same angu- 
lar velocity |^8|. Any destabiUzing, relativistic effect in 
binaries therefore has to be smaller. 

The collapse of binary neutron stars to black holes 
prior to merger reported by WMM could, in principle, 
be caused either by a decrease of the maximum allowed 
mass, or by a dynamical instability. As we have shown, 
the maximum allowed mass, within our assumptions and 
approximations, increases, which rules out the first pos- 
sibility. Note, however, that we are only constructing 
quasiequilibrium configurations, which may not be dy- 
namically stable. In |Q we show that all inspiraling 
binary neutron stars are secularly stable against radial 
collapse to black holes all the way down to the innermost 
stable circular orbit (ISCO). While this does not com- 
pletely rule out the existence of a dynamical instability, 
we note that in Newtonian binaries, dynamical instabili- 
ties always occur later along equilibrium sequences than 
secular instabilities Q,^. The same result has been 
shown for single, rotating relativistic stars [Q. Recently, 
Thorne |^ has argued analytically that tidal fields sta- 
bilize systems and that stars which are stable in isolation 
are stable with respect to both secular and dynamical 
modes in binary configurations. 

Fig. 2 demonstrates that at fixed rest-mass, the central 
density decreases as the stars approach each other and 
get tidally deformed. This effect, as well as the increase 
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FIG. 4. Blow-up of two curves in Fig. 2: binding en- 
ergy (solid) and angular momentum (dashed) as a function 
of the angular velocity for a binary with Mo = 0.169 and 
{M/R)a^ = 0.175. 



of the maximum allowed mass, is consistent with post- 
Newtonian predictions p3| , p^ , p^ . 

Next we construct sequences of constant rest-mass 
Mo, which up to the ISCO approximate evolutionary se- 
quences. As discussed in Section ||, we maintain coro- 
tation, whereas in reality it is more likely that circula- 
tion will be conserved. Nevertheless, our sequences are 
the first sequences of inspiraling binaries in full general 
relativity. Moreover, post-Newtonian sequences of con- 
stant circulation are not vastly different from corotat- 
ing sequences [15|. In Fig. 3 we plot the binding en- 
ergy {M — MooJJMo and the angular momentum J as a 
function of separation for several different rest-masses. 
Since the separation is not an invariant quantity, we 
have parametrized the sequence by the nondimension- 
alized angular velocity Mofl (= MqO.). Our curves do 
not connect to Mq^ = 0, corresponding to infinite sepa- 
ration, since we can numerically resolve only fairly close 
models. 

In the top half of Fig. 3 we show plots for sequences 
for several different, increasingly relativistic rest-masses. 
The curves are labeled by the compaction {M/R)oq that 
the stars would have in isolation at infinity. We have 
plotted graphs for {M/R)oo between 0.05 and 0.2 in in- 
crements of 0.025. In the lower half of Fig. 3 we show 
corresponding plots of J. According to Eq. ( |67| ) the min- 
ima in both curves must agree, which they do within our 
numerical accuracy. In Fig. 4 we show a blow-up of the 
two curves for stars with {M/R)ao — 0.175. 

For infinitely separated stars, both the binding energy 
and the angular velocity vanish. As the stars approach 
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Mo 


Moo 




MoQ,isco 


(Jtot/A^tot) /SCO 


0.059 


0.058 


0.05 


0.003 


1.69 


0.087 


0.084 


0.075 


0.0065 


1.37 


0.112 


0.106 


0.1 


0.01 


1.22 


0.134 


0.126 


0.125 


0.015 


1.12 


0.153 


0.142 


0.15 


0.02 


1.05 


0.169 


0.155 


0.175 


0.025 


1.00 


0.178 


0.162 


0.2 


0.03 


0.97 



TABLE II. Numerical values for sequences of constant 
rest-mass Mo and polytropic index n = 1. We tabulate the 
total energy Moo and compaction {M/R)oo each star would 
have in isolation as well as the angular velocity Mofl and the 
angular momentum Jtot/M^t the ISCO. The maximum 
rest-mass in isolation is M™'"' = 0.180. 



each other, the angular velocity increases while the bind- 
ing energy decreases. This effect is essentially Newtonian 
and is even evident for two Newtonian point masses. As 
the stars approach, however, finite size effects eventu- 
ally play an important role. The energy associated with 
the rotation of the individual stars adds to the (nega- 
tive) binding energy, and therefore reduces it. For stiff 
enough equations of state, for which the moment of in- 
ertia and hence the rotat iona l energy of the individual 
stars is large (see Section VB ), the binding energy goes 
through a minimum and then increases again prior to 
contact. The location of the minimum marks the onset 
of a secular instability, beyond which the binary can no 
longer maintain corotation. It is expected that the dy- 
namical instability defining the ISCO occurs after, but 
close to the onset of the secular instability 0,|l^ . In the 
following we will refer to the location of the minimum as 
the ISCO. 

The upper labels give the orbital frequency in Hz for 
stars of rest- mass 1.5Mq. The corresponding gravita- 
tional wave frequency is larger by a factor of 2 for the 
dominant quadrupole mode. For small values of the com- 
paction we find ISCO frequencies comparable to those 
reported by WMM. However, for larger compaction and 
more relativistic configurations we find frequencies very 
similar to what is found from post-Newtonian calcula- 
tions Hi. 

We summarize our results in Table 2, where we also in- 
clude the dimensionless angular momentum Jtot/Aftot ^ 
J/2M^ at the ISCO. For small rest-masses, this value 
is larger than unity, in agreement with WMM. For high 
enough rest-masses, however, it drops below unity, so 
that the two stars could plunge and form a Kerr black 
hole without having to lose additional angular momen- 
tum. 
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FIG. 5. Rest-Mass Mo of a n = 1.5 polytrope versus maxi- 
mum density pc for separations za = 0.3 (bottom solid line), 
0.2, 0.1 and 0.0 (top line). The dashed line is the Oppen- 
heimer-Volkoff result. 



B. Sequences for n = 1.5 and n = 2.0 

In this Section we will present results for polytropic 
indices of n = 1.5 and 2.0, representing softer equations 
of state. Except for the absence of an ISCO prior to con- 
tact (see below) all results are qualitatively very similar 
to those for n = 1. In particular, we consistently find a 
decrease of the maximum density as the stars approach, 
and an increase in the maximum allowed mass. The rel- 
ative size of these effects differs for three basic reasons: 
First, for softer equations of state, the maximum mass of 
a star occurs at a smaller value of the compaction M/R, 
and hence relativistic effects play a smaller role. Second, 
for softer equations of state these stars are more centrally 
condensed. We therefore expect tidal fields to play a less 
important role for the stability of these stars in close bi- 
naries. While it is easier to deform their surface, the 
bulk of the matter is very concentrated at the core of the 
stars and well shielded from the tidal field of the compan- 
ion. Third, for softer equations of state, the stars have a 
smaller orbital frequency even at very small separations, 
so that the effects of rotation are smaller. Accordingly 
we find that the maximum allowed mass still increases 
with decreasing separation, but the effect is smaller than 
for n = 1. 

More centrally condensed stars have a smaller moment 
of inertia, and hence the rotational kinetic energy asso- 
ciated with the spin of the stars is smaller than for less 
centrally condensed stars. Therefore a turning point in 
the binding energy curve can only be expected for stars 
with a stiff enough equation of state. This effect has been 
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FIG. 6. Binding energy of n = 1.5 polytropes as a func- 
tion of the angular velocity for different rest-masses. The 
curves are labeled by the compaction {M/R)oo of the stars 
in isolation. The maximum compaction for a stable, isolated, 
nonrotating n = 1.5 polytrope is 0.136. 

discussed by several authors in the context of Newtonian 
theory [ ^7|jl7|j4^ . We did not see a turning point for 
n > 1.5, in agreement with For these polytropic 

indices there is no ISCO prior to contact, and we expect 
the orbits to be stable until the stars touch and form a 
contact binary. This is the first construction of a contact 
binary in full general relativity. Proving the existence of 
a contact binary neutron star (by, e.g., the signature of 
its gravitational waveform) would indicate that the equa- 
tion of state of nuclear matter is rather soft. We do not 
expect this to be the case 

In Fig. 5 we plot the rest-mass versus the central den- 
sity for several different separations for n — 1.5. Quali- 
tatively the result is very similar to Fig. 2 for n = 1: For 
all separations the curve differs from the OV result by 
less than 1%. For decreasing separation we find a small 
increase in the allowed mass that a given density can 
support. In particular, the maximum quasiequilibrium 
rest-mass increases by roughly 1.2% from M = 0.275 for 
ZA = 0.3 to 0.278 for stars in contact. For n = 1 the cor- 
responding increase is about 2 %. The maximum density 
decreases as the stars approach and get tidally deformed. 

In Fig. 6 we plot the binding energy of n = 1.5 poly- 
tropes as a function of the angular velocity. We show 
results for several different rest-masses and label them 
by the compaction (M/i?) oo for the same stars in isola- 
tion. In contrast to the results for n = 1, these curves no 
longer show a turning point. This implies that the stars 
are secularly stable all the way to touching. 

In Figs. 7 and 8 we show the corresponding results for 
n = 2 polytropes. Again, in Fig. 7 we show the rest-mass 
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FIG. 7. Same as Fig. 5 for a n = 2.0 polytrope. 
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FIG. 8. Same as Fig. 6 for a n = 2.0 polytrope. The maxi- 
mum compaction {M/R)ao for a stable, isolated, nonrotating 
n — 1.5 polytrope is 0.075. 
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versus central density. The maximum quasiequilibrium 
rest-mass increases from M = 0.523 for za = 0.3 to 0.528 
for touching stars. This relative increase of roughly 1 % is 
smaller than even for n = 1.5. As expected, the binding 
energies in Fig. 8 do not show a turning point, so that 
the binaries are secularly stable all the way to touching. 

VI. SUMMARY AND CONCLUSIONS 



NASA Grant NAG 5-3420 at Ilhnois, NSF Grant PHY 
94-08378 at Cornell, and by the NSF Binary Black Hole 
Grand Challenge Grant Nos. NSF PHY 93-18152/ ASC 
93-18152 (ARPA supplemented). Computations were 
performed at the Cornell Center for Theory and Simula- 
tion in Science and Engineering and the National Center 
for Supercomputing Applications, University of Illinois 
at Urbana-Champaign. 



We report on the first fully relativistic calculation of 
binary neutron stars in quasiequilibrium. We previously 
presented some of our preliminary results in |2^ ; here 
we describe in detail all our assumptions and approx- 
imations, equations and numerical algorithm, as well 
as results for different polytropic indices. We integrate 
a subset of Einstein's equations, coupled to the equa- 
tion of hydrostatic equilibrium, to solve the initial value 
problem for binaries. We construct models of coro- 
tating binary neutron stars in close circular orbit, in- 
cluding relativistic models of contact binaries. We also 
construct sequences of constant rest-mass configurations 
parametrized by their separation and orbital angular fre- 
quency. 

We find that the maximum density of the stars de- 
creases as the stars approach and get tidally deformed. 
Simultaneously, the mass that a given maximum density 
can support increases as the stars approach. In particu- 
lar, we find that the maximum allowed mass of neutron 
stars in quasiequilibrium binaries increases with decreas- 
ing separation. These effects are larger for smaller poly- 
tropic index (and hence stiffer equation of state). 

Searching for turning points of the binding energy of 
constant rest mass sequences, we locate, for small enough 
polytropic index, the ISCO. As in the case of Newtonian 
configurations, an ISCO exists only for indices n ~ 1.5; 
for softer equations of state, contact is reached prior to 
the onset of orbital instability. 

In ||2^ we presented a more careful analysis of the ra- 
dial stability of relativistic binary neutron stars against 
collapse. We showed that all inspiraling binary neutron 
stars are secularly stable against radial collapse to black 
holes all the way to the ISCO (or contact, if, for large 
enough n, no ISCO is encountered). We do not find any 
evidence for a destabilization of neutron stars in close 
binary orbits. 
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APPENDIX A: NUMERICAL RESULTS FOR 
SELECTED SEQUENCES 

In the following we tabulate numerical values for se- 
lected sequences. For a given polytropic index n and the 
rest-mass (baryon mass) Mq of one star (or equivalently 
its compaction in isolation (M/i?)oo), we list the relative 
separation za — TA/rs, the maximal density parameter 
^max^ the mass M the angular momentum J, the (or- 
bital) frequency f2, and the locations fA, rs and fc- We 
have "barred" these quantities as a reminder that they 
are dimensionless coordinate values. Recall that p, po 
and P may be obtained from q via Eqs. d?]) - (||). 



„ = 1, Mo = 0.0595, (M/R)oo = 0.05 





max 

y 


M 


J 


n 


rA 


rs 


rc 


0.00 


0.0275 


0.057806 


0.01095 


0.061 


0.000 


1.529 


2.773 


0.10 


0.0278 


0.057806 


0.01094 


0.057 


0.281 


1.594 


2.810 


0.15 


0.0281 


0.057809 


0.01098 


0.053 


0.430 


1.677 


2.868 


0.20 


0.0284 


0.057815 


0.01109 


0.048 


0.591 


1.791 


2.959 


0.25 


0.0286 


0.057825 


0.01129 


0.042 


0.771 


1.940 


3.087 


0.30 


0.0288 


0.057836 


0.01155 


0.037 


0.975 


2.118 


3.251 








TABLE in. 










n 


= 1, Mo = 


= 0.1118, (M/R) 


CO =0.1 




ZA 


^max 


M 


J 


Cl 


rA 


rs 


fc 


0.00 


0.0658 


0.105511 


0.02715 


0.101 


0.000 


1.289 


2.353 


0.10 


0.0667 


0.105502 


0.02707 


0.094 


0.238 


1.346 


2.384 


0.15 


0.0676 


0.105509 


0.02710 


0.087 


0.365 


1.418 


2.433 


0.20 


0.0685 


0.105521 


0.02729 


0.079 


0.502 


1.516 


2.511 


0.25 


0.0693 


0.105558 


0.02766 


0.070 


0.655 


1.644 


2.621 


0.30 


0.0698 


0.105593 


0.02818 


0.062 


0.828 


1.797 


2.763 








TABLE IV. 










n — 


= 1, Mo = 


0.1341, {M/R)oo = 0.125 




ZA 


max 


M 


J 




rA 


rs 


fc 



0.00 0.0912 0.124786 0.03496 0.122 0.000 1.172 2.148 

0.10 0.0926 0.124785 0.03482 0.114 0.217 1.225 2.175 

0.15 0.0940 0.124787 0.03482 0.106 0.332 1.291 2.219 

0.20 0.0954 0.124819 0.03500 0.096 0.458 1.381 2.290 

0.25 0.0967 0.124849 0.03538 0.086 0.597 1.498 2.390 

0.30 0.0976 0.124898 0.03596 0.076 0.756 1.639 2.520 



TABLE V. 
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n = 1, Mo = 0.1534, {M/R)oo = 0.15 



ZA 




M 


J 




TA 


TB 


re 


0.00 


0.1235 


0.140851 


0.04188 


0.146 


0.000 


1.056 


1.943 


0.10 


0.1256 


0.140842 


0.04167 


0.137 


0.196 


1.104 


1.967 


0.15 


0.1280 


0.140846 


0.04162 


0.127 


0.300 


1.163 


2.005 


0.20 


0.1303 


0.140859 


0.04174 


0.116 


0.413 


1.244 


2.067 


0.25 


0.1325 


1 4090S 


0.04210 


0.104 


0.539 


1.350 


2.156 


0.30 


0.1341 


0.140971 


0.04268 


0.092 


0.682 


1.477 


2.273 








TABLE VI. 










n = 


1, Mq = 


0.1685, (M/i?)oo = 0.175 




ZA 




M 


J 


fi 


VA 


rB 


fc 


0.00 


0.1647 


0.152893 


0.04719 


0.173 


0.000 


0.944 


1.743 


0.10 


0.1683 


0.152883 


0.04691 


0.163 


0.176 


0.987 


1.762 


0.15 


0.1726 


0.152875 


0.04677 


0.152 


0.268 


1.038 


1.792 


0.20 


0.1769 


0.152893 


0.04681 


0.139 


0.368 


1.108 


1.844 


0.25 


0.1811 


0.152936 


0.04708 


0.125 


0.480 


1.201 


1.920 


0.:i() 


O.LSll 


().1.")2997 


l).()175« 


0.111 


().()()() 


l.:!i2 


2.022 








TABLE VII. 










n = 


= 1, Mq = 


= 0.1781, (M/i?) 


oo =0.2 




ZA 




M 


J 


fi 


fA 


fs 


fc 


0.00 


0.2164 


0.160183 


0.05024 


0.202 


0.000 


0.841 


1.560 


0.10 


0.2228 


0.160174 


0.04989 


0.191 


0.157 


0.877 


1.572 


0.15 


0.2327 


0.160137 


0.04963 


0.180 


0.238 


0.917 


1.587 


0.20 


0.2450 


0.160130 


0.04948 


0.168 


0.323 


0.970 


1.616 


0.25 


25Qn 


0.160145 


0.04953 


0.154 


0.415 


1.038 


1.662 


0.30 


0.2741 


0.160189 


0.04975 


0.139 


0.517 


1.119 


1.725 






TABLE VIII. 










n = 


1.5, Mo 


= 0.241, 


(M/ii), 


oo = 0.85 




ZA 


^max 


M 


J 


fi 


fA 


fs 


fc 


0.00 


0.0626 


0.231583 


0.13408 


0.035 


0.000 


3.409 


6.227 


0.10 


0.0633 


0.231623 


0.13471 


0.032 


0.631 


3.569 


6.318 


0.20 


0.0650 


0.231708 


0.13738 


0.027 


1.328 


4.014 


6.642 


0.30 


0.0665 


0.231853 


0.14341 


0.021 


2.184 


4.740 


7.281 



TABLE IX. 



n = 1.5, Mo = 0.258, {M/R)^ = 0.1 



ZA 


^max 


M 


J 




TA 


tb 


rc 


0.00 


0.0794 


0.246547 


0.14275 


0.042 


0.000 


3.038 


5.561 


0.10 


0.0802 


0.246600 


0.14346 


0.039 


0.564 


3.183 


5.643 


0.20 


0.0831 


0.246688 


0.14574 


0.033 


1.180 


3.564 


5.902 


0.30 


0.0855 


0.246887 


0.15166 


0.026 


1.938 


4.204 


6.460 



TABLE X. 





n = 


1.5, Mo = 


0.2745, (M/i?) 


oo = 0.125 




ZA 


^max 


M 


J fi 


fA 


fB 


rc 


0.00 


0.1119 


0.260578 


0.14820 0.055 


0.000 


2.549 


4.680 


0.10 


0.1141 


0.260614 


0.14830 0.051 


0.472 


2.658 


4.722 


0.20 


0.1237 


0.260665 


0.14908 0.045 


0.964 


2.909 


4.824 


0.30 


0.1380 


0.260810 


0.15227 0.038 


1.518 


3.291 


5.061 



TABLE XI. 





n 


= 2, Mo : 


= 0.495, (M/i?) 


oo = 0.05 




ZA 


^max 


M 


J n 


fA 


fB 


rc 


0.00 


0.0381 


0.48628 


0.7204 0.0073 


0.000 


12.89 


23.47 


0.10 


0.0383 


0.48635 


0.7281 0.0068 


2.382 


13.49 


23.82 


0.20 


0.0389 


0.48649 


0.7508 0.0057 


4.998 


15.11 


24.99 


0.30 


0.0395 


0.48672 


0.7924 0.0045 


8.192 


17.78 


27.30 



TABLE XII. 





n 


= 2, Mo ^ 


= 0.52, (M/ii)oo 


= 0.065 




ZA 


^max 


M 


J a 


fA 


fB 


rc 


0.00 


0.0493 


0.50929 


0.7190 0.0095 


0.000 


10.70 


19.52 


0.10 


0.0497 


0.50936 


0.7255 0.0089 


1.976 


11.18 


19.76 


0.20 


0.0525 


0.50950 


0.7409 0.0077 


4.069 


12.30 


20.34 


0.30 


0.0574 


0.50969 


0.7679 0.0065 


6.432 


13.95 


21.44 



TABLE XIII. 
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